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Gravitational lensing is a powerful astrophysical and cosmological probe and 
is particularly valuable at submillimeter wavelengths for the study of the sta- 
tistical and individual properties of dusty starforming galaxies. However the 
identification of gravitational lenses is often time-intensive, involving the sift- 
ing of large volumes of imaging or spectroscopic data to find few candidates. 
We used early data from the Herschel Astrophysical Terahertz Large Area 
Survey to demonstrate that wide-area submillimeter surveys can simply and 
easily detect strong gravitational lensing events, with close to 100% efficiency. 

When the light from a distant galaxy is deflected by a foreground mass — commonly a mas- 
sive elliptical galaxy or galaxy cluster or group — its angular size and brightness are increased, 
and multiple images of the same source may form. This phenomenon is commonly known as 
gravitational lensing (1) and can be exploited in the study of high-redshift galaxy structures 
down to scales difficult (if not impossible) to probe with the largest telescopes at present (2^) 
and to detect intrinsically faint objects. Surveys conducted at submillimeter wavelengths can 
particularly benefit from gravitational lensing because submillimeter telescopes have limited 
spatial resolution and consequently high source confusion, which makes it difficult to directly 
probe the populations responsible for the bulk of background submillimeter emission (5,6). In 
addition, galaxies detected in blank-field submillimeter surveys generally suffer severe dust ob- 
scuration and are therefore challenging to detect and study at optical and near infra-red (NIR) 
wavelengths. By alleviating the photon starvation, gravitational lensing facilitates follow-up ob- 
servations of galaxies obscured by dust and in particular the determination of their redshift (7). 
Previous submillimeter searches for highly magnified background galaxies have predominantly 
targeted galaxy cluster fields (8). In fact, a blind search for submillimeter lensing events re- 
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quires large area because of their rarity, and sub-arcseconds angular resolutions to reveal multi- 
ple images of the same background galaxies. Although the first requirement has recently been 
fulfilled, thanks to the advent of the South Pole Telescope (SPT) (9) and the Herschel Space 
Observatory (Herschel) (10), the second is still the prerogative of ground-based interferometric 
facilities, such as the Submillimeter Array (SMA) and the IRAM Plateau de Bure Interferom- 
eter (PdBI), which because of their small instantaneous field of view are aimed at follow-up 
observations rather than large-area survey campaigns. Nevertheless, several authors (11-14) 
have suggested that a simple selection in flux density, rather than surveys for multiply-imaged 
sources, can be used to easily and efficiently select samples of strongly gravitationally-lensed 
galaxies in wide-area submillimeter and millimeter surveys. The explanation for this lies in the 
steepness of the number counts (the number of galaxies at a given brightness) of dust-obscured 
star-forming galaxies, which are usually referred to as submillimeter galaxies (SMGs) (15). Be- 
cause of that, even a small number of highly-magnified SMGs can substantially affect the shape 
of the bright end of the submillimeter source counts enhancing the number of SMGs seen at 
bright flux densities than would be expected on the basis of our knowledge of the un-lensed 
SMG population (Fig.l). Furthermore, the frequency of lensing events is relatively high in 
the submillimeter (11) because SMGs are typically at high redshift (z >~ 1) (16), and this 
increases the probability that a SMG is in alignment with, and therefore lensed by, a fore- 
ground galaxy. Other important contributors to the bright tail of the submillimeter counts are 
low-redshift (z < 0.1) spiral and starburst galaxies (17) and higher redshift radio-bright Ac- 
tive Galactic Nuclei (AGNs) (18); however both of these are easily identified, and therefore 
removed, in relatively shallow optical and radio surveys. Therefore, flux-density limited sub- 
millimeter surveys could provide a sample of lens candidates from which contaminants can be 
readily removed, leaving a high fraction (close to 100 %) of gravitational lens systems (Fig. 1). 
Because this selection of lens candidates relies only on the properties of the background source 
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(its flux density), it can probe a wide range of lens properties (such as redshifts and masses) and 
thus provide a valuable sample for studying the elliptical properties of lensing galaxies (19) as 
well as investigating the detailed properties of the lensed SMGs. 

The submillimeter lens candidate selection at work. Although the approach presented above 

may be more efficient and vastly more time-effective than those exploited so far in the radio (20) 
or the optical (21, 22), at least several tens of square degrees (deg^) of the sky must be observed 
in the submillimeter to produce a statistically significant sample of strongly lensed objects 
and a minimal contamination from unlensed galaxies. This is because the surface density of 
lensed submillimeter galaxies is predicted to be lower than ~0.5 deg~^, for flux densities above 
100 mJy at 500 //m (Fig. 1). Submillimeter surveys conducted before the advent of Herschel 
were either limited to small areas of the sky (15, 23), or were severely affected by source con- 
fusion due to poor spatial resolution (24). Therefore no previous test of this selection method 
has been performed, although the SPT has recently mapped an area of more than 80 deg^ at 
millimetre wavelengths (9) and found an "excess" of sources that could be accounted for by a 
population of gravitationally-lensed objects. 

The Herschel Astrophysical Terahertz Large Area Survey (H- ATLAS) (25) represents the largest- 
area submillimeter survey being currently undertaken by Herschel. H- ATLAS uses the Spectral 
and Photometric Imaging REceiver (SPIRE) (26, 27) and the Photodetector Array Camera and 
Spectrometer (PACS) (28, 29) instruments and, when completed, will cover ~550 deg^ of the 
sky from 100 to 500 jira. H- ATLAS has been designed to observe areas of the sky with previ- 
ously existing multi-wavelength data: Galaxy Evolution Explorer (GALEX) ultra violet (UV) 
data, Sloan Digital Sky Survey (SDSS) optical imaging and spectroscopy, NIR data from the UK 
Infrared Telescope (UKIRT) Infrared Deep Sky Survey (UKIDSS) Large Area Survey (LAS), 
spectra from the Galaxy And Mass Assembly GAMA (30) project, radio imaging data from the 
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Faint Images of the Radio Sky at Twenty-cm (FIRST) survey and the NRAO Very Large Array 
Sky Survey (NVSS). The first 14.4 deg^ of the survey, centred on J2000 RA 09:05:30.0 DEC 
00:30:00.0 and covering ~3% of the total area, was observed in November 2009 as part of the 
Herschel Science Demonstration Phase (SDP). The results were a catalog of ~6600 sources (31) 
with a significance > 5 cr, in at least one SPIRE waveband, where the noise, a, includes both 
instrumental and source confusion noise and corresponds to ~7 to 9 mJy/beam. 
The Herschel/SPIRE 500 /im channel is favourable for selecting lens candidates, because the 
submillimeter source counts steepen at longer wavelengths (24, 32). We used theoretical pre- 
dictions {14) to calculate the optimal limiting flux density, above which it is straightforward 
to remove contaminants from the parent sample and maximize the number of strongly lensed 
high-redshift galaxies. The surface-density of un-lensed SMGs is predicted to reach zero by 
5*500 ~ 100 mJy {14) and these objects are only detectable above this threshold if gravitation- 
ally lensed by a foreground galaxy (Fig. 1). The H- ATLAS SDP catalog contains 11 sources 
with 500 lira flux density above 100 mJy. Ancillary data in the field revealed that six of these 
objects are contaminants; four are spiral galaxies with spectroscopic redshifts in the range of 
0.01 to 0.05 [see {33) for a detailed analysis of one of these sources], one is an extended Galac- 
tic star forming region, and one is a previously known radio— bright AGN {34). Although the 
number of these sources are few at bright flux densities, the measured surface densities are con- 
sistent with expectations (Fig. 1) {17, 18). Exclusion of these contaminants left the five objects 
that form our sample of lens candidates (table SI) {35), identified as ID9, IDl 1, ID17, ID81 and 
ID130. 

Unveiling the nature of the lens candidates. For gravitational lensing systems selected at 
submillimeter wavelengths, we would expect the lensing galaxy to be seen in optical and/or 
NIR images, in which the emission from the lens dominates over the higher redshift background 
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SMG. In line with these expectations, all of the lens candidates have a close counterpart in SDSS 
or UKIDSS images (or both). A likelihood ratio analysis (36) showed that the probability of a 
random association between these bright submillimeter sources and the close optical/NIR coun- 
terparts is less than a few per cent. Therefore the optical and submillimeter emissions must be 
physically related, either because they occur within the same object or because of the effects 
of gravitational lensing, boosting the flux of the background source and indirectly affecting the 
likelihood ratio calculations. The redshift measurements support the later scenario. Although 
the optical/NIR photometric/spectroscopic redshifts lie in the range of 2; ~ 0.3 to 0.9 (table 
1 and figs S3 and S4) (35), the redshifts estimated from the submillimeter/millimeter spectral 
energy distributions (SED) (table 2) [following the method described in (37, 38)], are distinctly 
different (Table 1). The lensed SMG photometric redshifts have been confirmed and made 
more precise through the spectroscopic detection, in these objects, of carbon monoxide (CO) 
rotational line emission which are tracers of molecular gas assocatiated to star forming environ- 
ments. Until recently, these kind of detections were difficult to achieve without prior knowledge 
of the source redshift, which required extensive optical/NIR/radio follow-up observations. Be- 
cause of the development of wide-bandwidth radio spectrometers capable of detecting CO lines 
over a wide range of redshifts, it is now possible for blind redshift measurements of SMGs to 
be taken without relying on optical or NIR spectroscopy (39, 40). ID81 was observed with the 
Z-Spec spectrometer (41, 42) on the California Institute of Technology Submillimeter Obser- 
vatory (CSO). The data revealed several CO lines redshifted into the frequency range of 187 
to 310 GHz; the strongest of these lines has been interpreted as the CO J=7— 6 line, with an 
estimated redshift of ^ = 3.04 (43). This represents the first blind redshift determination by 
means of Z-Spec. We followed up this observation with the PdBI and detected CO J=3— 2 
and CO J=5— 4 emission lines, redshifted io z = 3.042, confirming the Z-Spec measured red- 
shift (35). We also used the Zpectrometer instrument (44, 45) on the NRAO Robert C. Byrd 
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Green Bank Telescope (GBT) to obtain an independent confirmation of the redshift of ID81 
(table 1 and fig. SI) (35, 46) and to measure the redshift of ID130. We detected redshifted CO 
J=l— emission sA. z — 2.625 in the spectrum of ID130 (fig. SI) {35, 46)]. This redshift was 
confirmed by the PdBI with the observation of CO J=3-2 and CO J=5-4 lines, yielding a redshift 
of z = 2.626 (35). The Z-Spec spectrometer observed the remaining three lens candidates (43) 
and detected CO lines at redshifts of z = 1.577 and z = 1.786 for ID9 (fig.S2) (35) and 
DDI 1, respectively, which are higher and inconsistent with the redshifts derived from the optical 
photometry /spectroscopy (Table 1). The Z-Spec CO measurements for ID 17 are indicative of 
two redshifts; one, z — 0.942, that is in agreement with the optical redshift and a higher one, 
z — 2.31, which is indicative of a more distant galaxy. 

To determine the morphological type of the foreground galaxies we obtained high resolution 
optical images for all five objects with the Keck telescope at g- and i-bands (35). ID9, IDU, 
ID81 and ID 130 all have optical profiles that are consistent with elliptical galaxies (figs S5 and 
S6 and table S4) (35). The interpretation of the results for ID 17 is complicated by the presence 
of two partially superimposed galaxies in the optical images (fig. S7) (35), neither exhibiting 
the disturbed morphology expected for lensed objects. This indicates that ID 17 may be a grav- 
itational lens system with two foreground lensing masses at similar redshifts (z ~ 0.8 to 0.9) 
— possibly a merging system — with some molecular gas responsible for the CO emission de- 
tected by Z-Spec at 2; ~ 0.9 and confirmed with optical spectroscopy (table 1). A fit to the 
UV/optical/NIR SEDs of ID9, IDll, ID81 and ID130 (47), using the models of (48), gives stel- 
lar masses in the range of 4 x 10^° to 15 x Mq (Table 2) and almost negligible present-day star 
formation, which is consistent with elliptical galaxies (fig. 2). 

For all five lens systems the background source appears to be undetected in the Keck g- and 
i-band images, despite the flux magnification due to lensing. After subtracting the best fit light 
profile from each lens we found no structure that could be associated with the background 
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source in the residual images (figs S5 and S6) (55). We derived S-cr upper limits from the 
residual maps (table S4) (35) and corresponding NIR limits from the UKIDSS images. These 
upper limits were used to fit the SEDs of the background sources assuming the models of (48), 
calibrated to reproduce the UV-to-infrared SEDs of local, purely star-forming ultra luminous 
infrared galaxies (ULIRG; lO^^ <Ljr/L0 <10^'^) (49) (fig. 2). A visual extinction (50) ofAy 
>2 is required to be consistent with the optical/NIR upper limits (fig. 2 and table 2), confirming 
severe dust obscuration in these galaxies along the line-of-sight. Our results indicate that these 
submillimeter bright gravitationally lensed galaxies would have been entirely missed by stan- 
dard optical methods of selection. 

We obtained observations at the SMA for ID81 and ID 130 at 880 /xm, with the aim of detect- 
ing the lensed morphology of the background galaxy (35). The SMA images reveal extended 
submillimeter emission distributed around the cores of the foreground elliptical galaxies, with 
multiple peaks (four main peaks in ID81 and two in ID 130), which is consistent with a lensing 
interpretation of these structures (Fig. 3). The position of these peaks can be used to directly 
constrain the Einstein radius — the radius of the circular region on the sky (the Einstein ring) 
into which an extended source would be lensed if a foreground galaxy were exactly along the 
line of sight of the observer to the source (for a perfectly circular lens). The Einstein radius is a 
measure of the projected mass of the lens, so it can be used to derive the total (dark plus lumi- 
nous) mass of the galaxy within the Einstein radius (table 2) (35). Another measure of the total 
mass of a lens is the line-of-sight stellar velocity dispersion, cr^. We have estimated from 
the local Faber— Jackson (FJ) relation (51) between ay and the rest-frame B-band luminosity 
for elliptical galaxies. Assuming passive stellar evolution for the lens galaxies, which is appro- 
priate for elliptical galaxies, we have extrapolated their rest-frame K-band luminosity to z = 
[using the evolutionary tracks of (52)], and then converted this to B-band luminosity using the 
B - K = 4.43 color relation from (53). The result was then applied to the FJ relation from (54). 
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Given a mass model for the lens (35), we can predict the Einstein radius of the galaxy from the 
value of (Tv expected from the FJ relation and compare it with that directly measured from the 
SMA images (Table 2). Although the value of the Einstein radius derived from the line-of-sight 
stellar velocity dispersion is affected by large uncertainties (as a result of the scatter in the FJ 
relation) it is consistent with the value measured in the SMA images for both 1D81 and ID 130. 
In order to test whether the properties of the lensing galaxies in our sample are consistent with 
those of other known lens ellipticals at similar redshift, we compared the V-band mass-to-light 
ratio of the lens galaxy for 1D81 and ID 130 (Table 2) (35) to those measured in the Sloan Lens 
Advanced Camera for Surveys (ACS) (fig. 4) (55), which cover a similar redshift range (z~0. 1 
to 0.3). The agreement with the average trend revealed by SLACS confirms that our lens se- 
lection method is not biased to lensing ellipticals with atypical luminosities. Moreover, the 
location of ID 130 in Fig. 4 indicates that our selection method can probe lower masses and 
lower luminosity lens galaxies than those sampled by SLACS, thus offering a wider range in 
lens properties to be investigated. 

The best fit SED to the submillimeter/millimeter photometry for each of the five background 
sources give infrared luminosities Lir >~3 x 10^^ L© (Table 2), which would classify these ob- 
jects as Hyper Luminous Infra-Red galaxies (HLIRGs; Ljr >1O^^L0). However, a correction 
for magnification because of lensing will reduce these values by a factor of 10 or greater. For ex- 
ample, assuming that the light distribution of the background source is described by a Gaussian 
profile with a full width at half maximum (FWHM) of 0.2 arcseconds [which is consistent with 
the physical extension of the background galaxy in (4)], the best-fit lens model (fig. S9) (35) 
predicts a total amplifications of ~19 and ~6 for ID81 and ID 130, respectively. Typical ampli- 
fications of 8 to 10 are also suggested by (14), therefore, it is more likely that these sources are 
ULIRGs. 

These results already provide constraints for models of the formation and evolution of massive 
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galaxies at high redshift. The fact that many (if not all) of the brightest SMGs detected in the H- 
ATLAS SDP field are amplified by lensing, implies that un-lensed z > 1 star-forming galaxies 
with flux densities more than 100 mJy at 500 //m are rare, with < 4.6 of them per 14.4 deg~^, 
at 99% probability (assuming Poisson statistics). This translates into a 0.32 deg"^ upper limit 
on the surface density of these sources. The same limit should translate to the abundance of 
HLIRGs with Ljr >5x 10^^ Lo^t z < A, because they would also have 500-//m flux densities 
above 100 mJy, which has possible implications for the role of feedback during the formation 
of the most massive galaxies in the universe. By extrapolating our SDP findings to the full 
H- ATLAS field, we predict a total sample of more than 100 bright lensed sources, with which 
we can further improve this constraint. 
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Tables 



Table 1: Photometric and spectroscopic redshifts of the five lens candidates. Spectroscopic 
redshifts were derived from optical lines for the lens [^^sp^c^] and from CO lines for the back- 
ground source [^^pe^^]- Photometric redshifts are based on UV/optical/NIR photometry for the 
lens [4r*^] and H- ATLAS plus SMA and Max-Planck Millimeter Bolometer (MAMBO) pho- 

. .j;..., J r (submillimeter/millimeter) ■ . . . ■ j u-r^. j 

tometry for the background source [Zpj^^^^ ' , usmg the photometric redshift code 

of (37, 38)]. The quoted errors on the redshifts correspond to a 68% confidence interval. 



^phot 



(sub— mm/millimeter) 
^phot 



SDPID 
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11 
17 
81 



0.679±0.057 

0.72±0.16 
0.77±0.13 
0.334±0.016 



0.7932±0.0012(2) 
0.9435±0.0009(2) 
0.2999±0.0002(3) 



1 iTOTS" 

1 Q+0-4 
^•^-0.3 

7 Q+0.2 



130 0.239±0.021 0.2201±0.002(6) 



2.6 



+0.4 
-0.2 



1.577±0.008(i^ 

1.786±0.005(i) 
0.942±0.004& 2.308±0.01l(i) 
3.037±0.010(^) 
3.042±0.00l(*^'(5) 
2.625±0.00l('') 
2.6260±0.0003(5) 



Datum is from CSO/Z-Spec (43) 
Datum is from the William Herschel Telescope (55) 
(3) Datum is from SDSS 

Datum is from GBT/Zpectrometer (46) 
Datum is from PdBI (35) 

Datum is from the Apache Point Observatory (35) 
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Table 2: Derived parameters for the five lens candidates. Estimated mass in stars (M*) and Star 
Formation Rate (SFR) of the foreground galaxy derived from the best-fit to the UV/optical/near- 

IR part of the SED; the Einstein radius measured from the SMA images (^e); mass within 
the Einstein radius (Me) estimated from ^e; line-of-sight stellar velocity dispersion (cr^'^) de- 
rived from the Faber-Jackson relation and the B-band luminosity produced by the best-fit to 
the UV/optical/NIR SED; Einstein radius (9^"^) calculated from a^"^; infrared luminosity of the 
background source (Lir), without correction for magnification, derived by fitting the submil- 
limeter/millimeter part of the SED and the upper limits at optical and NIR wavelengths (Fig. 2); 
and visual extinction parameter (Ay) inferred for the background galaxy. All the quoted errors 
correspond to a 68 per cent confidence interval. For ID 17 only the infrared luminosity and the 
extinction parameter of the background source are quoted because the lensing mass probably 
consists of two galaxies that can only be disentangled in the Keck images. The symbols M© 
and Lq denote the total mass and the total luminosity of the Sun, respectively, and correspond 
to Mq = 1.99 X 10^° kg and Lq = 3.839 x 10^^ erg s~^. Dashes indicate lack of constraints. 
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Figure 1: Selection of gravitational lenses at submillimeter wavelengths. The 500 /im source 
counts consist of three different populations {14): high-redshift SMGs; lower redshift late type 
(starburst plus normal spiral) galaxies; and radio sources powered by active galactic nuclei. 
Strongly lensed SMGs dominate over unlensed SMGs at very bright fluxes where the count 
of un-lensed SMGs falls off dramatically (yellow shaded region). The data points are from 
H-ATLAS {32). 
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Figure 2: Spectra of the gravitational lens candidates. The UV, optical and NIR data points (blue 
dots) are from GALEX, SDSS and UKIDSS LAS, respectively. The submillimeter/millimeter 
data points (red dots) are from PACS/Herschel, SPIRE/Herschel, SMA and Max-Planck Mil- 
limeter Bolometer (MAMBO)/IRAM. Upper limits at PACS/Herschel wavelengths are shown 
at 3 a. ID 130 lies outside the region covered by PACS. The photometric data were fitted using 
SED models from (48). The background source, responsible for the submillimeter emission, is 
a heavily dust obscured star-forming galaxy (red solid curve), whereas the lens galaxy, which is 
responsible for the UV/optical and NIR part of the spectrum, is caracterized by passive stellar 
evolution. 
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Figure 3: Submillimeter and optical follow-up imaging of ID81 and ID 130. The SMA images 
of ID81 and ID 130 are shown in the top panels, centered on the optical counterpart, and were 
obtained by combining the visibility data from very-extended, compact and sub-compact con- 
figuration observations. The Keck i-band image of ID81 and ID130 are shown in the bottom 
panels with the SMA contours superimposed (in red). The contours are in steps of -2cr, 2a, 
4a, 6a, Sa, lOa... , with a = 0.6mJy/beam. The SMA synthesized beam is shown in the 
bottom-left comer. 
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Figure 4: Relationship between mass and luminosity for the lensing galaxy in ID81 and ID 130. 
The rest-frame band luminosity was derived from the best-fit SED to the UV/optical/NIR 
photometric data; the mass within the Einstein radius is that measured directly from the SMA 
images. The light versus mass relation inferred for ID81 and ID130 (red dots) is consistent with 
that observed for the SLACS lenses [black dots, from (55) assuming an uncertainty of 0.025dex 
in their mass estimates]. 
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Supporting Online Material 



MAMBO observations 

The five H-ATLAS/SDP lens candidates, i.e. ID9, IDll, ID17, ID81 and ID130, were ob- 
served on 2010 March 9 and 10, with MAMBO, at the IRAM 30 meter telescope on Pico 
Veleta, in Diretcor's Discretionary Time (DDT). The MAMBO array consists of 1 17 bolometer 
elements and operates at a central frequency of 250 GHz, corresponding to 1.2 mm. The beam 
size (FWHM fa 1 1 arcsec) of MAMBO ensures that the true dust emission at 1.2 mm is obtained 
if the source is not more extended than a few arcseconds. Each science target was observed in 
the photometric mode ("on-off") of MAMBO which is based on the chop-nod technique and 
placing the target on a reference bolometer element (on-target channel). The total observing 
time was 1.5 hours. Data were reduced using MOPSIC, and the current version of MOPSI 
(Zylka 1998, The MOPSI Cookbook (Bonn: MPIfR)). 

Keck observations 

The imaging observations were acquired on 10 March 2010 using the dual-arm Low Resolution 
Imaging Spectrometer [LRIS; (i, 2)] on the 10-m Keck I telescope. Each target field received 
simultaneous 110 x 3 s integrations with the s'-filter and 60 x 3s integrations with the i-filter 
using the blue and red arms of LRIS. A ~ 20" dither pattern was employed to generate on-sky 
flat field frames when incorporating all five fields. In addition, 1 s integrations were acquired in 
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the g- and i-filters for photometric caUbration of bright stars in each field. The seeing FWHM 
for the science exposures was ~ 0.8 — 0.9". The data were reduced using IDL routines and 
combined and analysed using standard IRAF tasks. 

SMA observations 

Observations of two H-ATLAS/SDP sources, ID81 and ID130, were obtained at 880 /xm us- 
ing the only current high resolution submillimeter facility in the world, the SMA. The SMA is 
an interferometer located near the summit of Mauna Kea, Hawaiim and consists of eight 6 m di- 
ameter radio telescopes. The two H- ATLAS sources were observed in Director's Discretionary 
Time (from February to May 2010) in three separate configurations, with baselines spanning a 
spatial range from 6 to 509 meters, over a total of 4 observing periods (Table S3). 
Target observations from each period were interspersed with observations of calibration sources, 
quasars J0909+013 and J0825+031 (phase) and J0730-116 and J0854+201 (amplitude). The 
phase calibration targets were typically observed every 7 to 15 minutes, depending upon con- 
figuration (a faster cycle was used for the larger configurations). 

Calibration of the complex visibility data was performed within the SMA's MIR package, a 
suite of IDL-based routines designed for use with SMA data. The initial opacity correction 
was obtained through application of system temperature to the raw visibility data, a standard 
practice. Further complex gain corrections, to remove both atmospheric and instrumental ampli- 
tude and phase variations, were measured using the calibration quasars, which appear as point 
sources to the interferometer. Calibration of the absolute flux density scale was performed us- 
ing measurements of Titan, whose continuum and line structure is known to within about 5% at 
submillimeter/millimeter wavelengths . 
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The resulting calibrated visibility data for each source were combined and imaged within the 
NRAO Astronomical Image Processing System (AIPS). Photometry obtained for the SMA im- 
ages along with those from PACS, SPIRE, MAMBO are given in table SI. 

Plateau de Bure Interferometer observations 

The H-ATLAS/SDP sources ID81 and ID 130 were observed in the CO J=3-2 and CO J=5-4 
lines with the IRAM Plateau de Bure Interferometer (i). Both sources were observed in excel- 
lent atmospheric conditions and with the full sensitivity of the six-element array. The observing 
frequencies were based on the redshifts provided by the CSO/Z-spec spectrometer. The receiver 
bandwidth was adjusted for maximum sensitivity and the observing frequencies centered in the 
1 GHz baseband of the narrowband correlator. Observations of ID81 were made on March 22, 
2010 for an effective integration time of 22 min and 14 min, respectively, for the CO J=3-2 and 
J=5-4 lines. The RF calibration was measured on 3C84, and amplitude and phase calibrations 
were made on 0823-1-033. The J=3-2 and J=5-4 transitions in ID130 were observed on March 
26 and April 16, 2010, respectively, for an effective integration time of 74 min and 32 min. 
The RF calibration was measured on 3C273, and amplitude and phase calibrations were made 
on 0906-1-015. The absolute flux calibration scales for ID81 and ID 130 were established using 
as primary calibrator MWC349. Data reduction and calibration were made using the GILDAS 
software package in the standard antenna based mode. 

Optical spectroscopic observations 
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Optical spectroscopic observations of ID 11 and ID 17 were made using the ISIS double-arm 
spectrograph on the 4.2-m William Herschel Telescope (WHT). The R158B and R158R grat- 
ings were used to provide wavelength coverage across the entire optical spectrum, split by 
a dichroic at ~5300 A. Four 900-second exposures were taken of each source in a standard 
'ABBA' pattern, nodding the telescope along the slit by 10 arcseconds between the first and 
second exposures, and back to the original position between the third and fourth integrations. 
This allowed initial sky subtraction to be performed by simply subtracting the 'A' frames from 
the 'B' frames. Additional sky subtraction was performed by subtracting the median value of 
each row, and then the positive and negative beams were aligned and coadded. Wavelength 
calibration was performed using observations of arc lamps taken with the same set-up. A one- 
dimensional spectrum was then optimally extracted. The spectra were taken through thin cloud 
and therefore no attempt has been made to flux-calibrate them. There was very little signal 
in the blue arms and so only the red-arm spectra are presented here. The redshifts of the two 
sources were determined by cross-correlation with template spectra. All reduction steps were 
undertaken using the IRAF package. The resulting spectra are shown in Fig. S3. The spec- 
trum of ID 1 1 reveals absorption lines associated to singly ionised calcium Ca H+K (rest-frame 
wavelengths: 3968.5 A for H-line and 3933.7 A for K-line) and the 4000 A break feature (rest- 
frame wavelength 4000 A) at z — 0.793, while the spectrum of ID 17 shows the emission from 
oxygen doublet [On]3727 (rest-frame wavelengths 3726-3729 A) and the 4000 A break feature 
at z = 0.944. In both spectra the absorption feature observed at ~ 7600 A is due to the Earth's 
atmosphere. 

A 30-minute exposure of ID130 was taken on May 15, 2010, with the Apache Point Observa- 
tory's 3.5-meter telescope and the DIS [Dual Imaging Spectrograph, (4)] long-slit spectrograph 
through medium clouds at an average airmass of 1.5. The spectrum was processed by sub- 
tracting the detector bias, dividing by a flat-field frame to correct for variable pixel response. 
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perfomiing distortion correction to align the spectrum in the wavelength and spatial directions, 
subtracting the sky flux determined from parts of the slit containing no sources, and applying 
a wavelength calibration by reference to emission lines from a Helium-Neon- Argon calibration 
lamp. Two emission lines in the spectrum (Fig. S4) were identified as [O n]3727 and [Ne 
III] 3 869 (rest-frame wavelengths 3869 A) from the ratio of their observed wavelengths. From 
the ratio of their observed to emitted wavelengths the redshift of the galaxy was determined to 
be z =0.2201±0.002. 
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Modelling with GALFIT 



GALFIT (5) is a publicly available two-dimensional non-linear fitting algorithm, which allows 
galaxy images to be modelled with one or multiple analytical light profiles. Each profile is 
constrained by a function and a set of parameters. GALFIT convolves the profiles with a user 
supplied point spread function, in this case empirical point spread functions constructed using 
nearby stars, and then performs a least-squares minimisation. No hard or soft constraints were 
applied to the fitting parameters to avoid any prior on the galaxy morphological type. For ID9 
and ID 11 single Sersic profiles resulted in a reduced close to 1.0 (see table S4 for the best 
fit parameters). ID17 was fitted with two Sersic component, assuming two lensing galaxies. 
The resulting Sersic indices were both less than 1 (see table S4). For ID81 and ID 130 two 
components were necessary to achieve a satisfactory fit, with a clean residual. The best fits 
were obtained using a combination of a compact elliptical Sersic core plus an exponential disk. 
No detectable background structure was revealed after subtracting the models, which shows the 
background galaxy is below the optical detection limit. Postage stamp images of ID9, IDll, 
ID17, ID81 and ID130 are shown in FigsS5 and S6, together with the corresponding best-fit 
models and residuals, while Fig. S7 shows the individual GALFIT components for IDS 1 and 
ID130. 

To derive photometric upper limits, we performed random aperture photometry on the i- and 
g-band Keck maps, using a 1.5 arcsecond radius. This radius was chosen to correspond with the 
structure visible in the SMA images for ID81 and ID 130, which extends to regions with radii 
of approximately 1 — 1.5 arcseconds. The resulting flux distributions were fitted with Gaussians 
and the 3— cr upper limits are presented in Table S4. 
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Mass estimate from lensing 



The Einstein radius of a strong galaxy-galaxy gravitational lens system can be measured from 
the configuration of multiple lensed images by averaging the distances of the images from the 
center of the lensing galaxy. For two of the H-ATLAS/SDP lens candidates, ID81 and ID 130, 
the positions of the lensed images are constrained by high-resolution SMA follow-up imaging. 
The lensed images of the background sources appear as peaks in the SMA signal-to-noise ratio 
map. Here we have selected those peaks with signal-to-noise ratio above eight, which pro- 
vided positions for four images in ID81 and two images in ID 130. The error on the Einstein 
radius is estimated by taking into account the uncertainties on the position of the individual 
peaks. For a point source the rms error on its position is A/2cr/SNR (assuming no systematic 
astrometry errors and uncorrelated Gaussian noise), where a is the Gaussian rms width of the 
instrument beam ( = FWTIM/2-\/21n2), while SNR is the signal-to-noise ratio at the source posi- 
tion (6, 7). The SMA synthesised beam (derived by combining observations in VEX, COM and 
SUB configurations) has size 0.81"x0.73" for ID81 and 0.78"x0.72" for ID130. Therefore, 
in estimating the relative positional uncertainty of the peaks, we have assumed FWHM=0.75" 
and FWHM=0.77" for ID81 and ID 130, respectively. The absolute positional uncertainty of 
the SMA images is estimated by referencing the data to nearby point-like sources (quasars) of 
known positions and is below 10 milli-arcseconds. 

Once the Einstein radius is known, the mass within the Einstein ring can be easily derived as- 
suming a Singular Isothermal Sphere (SIS) model (although the result is only little dependent 
on the model used) which is characterized by a projected surface density that falls off as 9~^, 
where 9 is the angular distance from the center of the mass distribution (8), 

Me = M(< ^e) = TTEerit^L (1) 
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and Ecrit is the critical surface density: 

E . - (2) 

In the equation above, c is the speed of Ught, G is the gravitational constant, Di, and Ds are 
the angular diameter distances to the lens and the source, respectively, while Dls is the angular 
diameter distance between the lens and the source. The error on the mass is obtained by prop- 
agating the errors on the Einstein radius and on the spectroscopic/photometric redshifts used to 
derived the angular diameter distances. The estimated values of 9e and Me are listed in Table 2. 
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Gravitational lensing modeling 



A detailed analysis of the lensed structure revealed by the SMA images is beyond the scope 
of this paper and is deferred to a forthcoming publication. However, in order to prove that such 
a structure is consistent with a lensing event, we have used the publicly available LENSMODEL 
software (9) to reproduce the positions of the peaks in the SMA maps. We have assumed a Sin- 
gular Isothermal Ellipsoid (SIE) model (8) for the mass distribution of the lens galaxy. The SIE 
model consists of concentric and aligned elliptical isodensity contours with axis ratio q. The 
circular limit is the SIS model and corresponds to g = 1. The results are shown in Fig. S9. We 
have further assumed that the centroid of the mass model coincides with that of the light distri- 
bution of the lensing galaxy. The best-fit model for ID 130 has ellipticity e = 0.16 and position 
angle (measured East of North) of ^ = +75 deg, consistent with the results found for the optical 
light-distribution that is dominated by the more compact Sersic profile (Table S4 and Fig. S8). 
For IDS 1, the mass distribution has ellipticity e = 0.24 and position angle 6 — —3 deg, which is 
not consistent with that measured for the luminous component (Table S4 and Fig. S8). Besides, 
the position of the peaks is not well reproduced by the model. This may hint at the effect of an 
external shear (which we did not include) due to a nearby cluster (photometrically detected 3.6 
arcminutes away), in the direction indicated by the arrow in Fig. S9. 

We have used the best-fit lens models to approximatively quantify the magnification experi- 
enced by a background source described by a Gaussian profile with a Full Width at Half Maxi- 
mum (FWHM) in the range 0.1-0.3". This extension is consistent with the physical size of the 
submillimeter galaxy studied by (10). The inferred magnification is ~18-31 for ID81 and ~5-7 
for ID 130. An example of lensed image, after convolution with the SMA point spread function, 
for the case FWHM=0.2" is shown in Fig. S9. 
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Figures 



Figure S 1: CO line detections in ID81 and ID 130. The figure shows the difference between 
the spectrum of ID81 and that of ID130, derived from Zpectrometer observations. The relative 
spectrum is normalized such that the peak line strength of ID81 is equal to 1. In both objects 
the peak is associated with the CO J=l-0 emission line (rest-frame frequency 1 15.27 GHz). 
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Figure S 2: CO line detections in ID9. The figure shows the spectrum derived from Z-Spec 
observations. The CO emission lines redshifted into the frequency range probed by Z-Spec cor- 
respond to transitions J=5-4 (rest-frame frequency 576.3 GHz) and J=6-5 (rest-frame frequency 
691.5 GHz). 
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Figure S 3: Optical spectra of ID 1 1 and ID 17 obatined with the WHT. 
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Figure S 4: Optical spectrum of ID 130 obtained with the APO telescope. The bottom panel of 
the figure shows the reduced 2-d spectrum in the region of the detected emission lines. The top 
panel shows the flux summed in a 5 -pixel wide (2 arcsecond) aperture centered on the object, 
with an arbitrary flux scale because the clouds made wavelength calibration impossible. 
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Figure S 5: Best-fit to the light distribution of the lens galaxy in the gravitational lens systems 
ID9, IDl 1 and ID 17. The postage stamp images show, from left to right, the keck i-band image, 
the best-fit light distribution model provided by GALFIT and the residual map obtained by sub- 
tracting the best-fit model from the observed light distribution. The map of the residuals show 
no evident structure, implying that the background source is particularly faint in the optical, 
despite the magnification due to lensing. 
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Figure S 6: Best-fit to the light distribution of the lens galaxy in the gravitational lens systems 
ID81 and ID 130. The postage stamp images show, from left to right, the keck i-band image, 
the best-fit light distribution model provided by GALFIT and the residual map obtained by 
subtracting the best-fit model from the observed light distribution. The SMA contours (in red) 
are overlaid on the optical images (in steps of 6a, 8cr, 10a, etc.) to highlight that there is 
no evident correspondence between the submillimeter and the optical emission in the residual 
maps. 



4 
3 
2 

c 
o 

Si 
< -1 

-2 
-3 
-A 



1 np — I 1 1 1 1 — 

A ' ID 81 Keck i 




J 



—I I I I I I L- 



1 1 1 1 1 1 1 

B GALFIT model 




—I I I I I I L- 



— 1 1 1 1 1 1 r- 



Residual . 




—I I I I I I i_ 



4 3 2 1 -1-2 -3 -4 3 2 1 0-1-2-3-4 3 2 1 0-1-2-3-4 

Arc Seconds 



4 r 1 1 1 1 1 1 1— 



— I 1 1 1 1 1 r- 



— 1 1 1 1 1 1 r- 



3 
2 

c 
o 

^ 
ui 

u 

< -1 

-2 
-3 
-A 



D 



ID130 Keck i 




GALFIT model 




—I I I I I I L- 



Residual 




—I I I I I I I— 



4 3 2 1 -1-2 -3 



-4 3 2 1 0-1-2-3-4 3 2 1 0-1-2-3-4 
Arc Seconds 



37 



Figure S 7: Decomposition of the best-fit models of the lens galaxies in ID17, ID81 and ID130. 
Top: ID 17 shows two partially superimposed components, indicative of two distinct lens galax- 
ies, each described by a relatively shallow Sersic profile. Middle: ID81 has one single lens 
galaxy whose light profile is reproduced by the sum of a Sersic profile and an exponential disk 
profile. Bottom: ID 130 is similar to ID81, with the light profile being described by the superpo- 
sition of a compact Sersic profile and an exponential disk profile. In both ID81 and ID 130, the 
SMA contours (in red) are overlaid on the optical images, in steps of 6cr, 8(T, lOcr, etc. 
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Figure S 8: Light profile of the lens galaxy in the gravitational lens systems ID81 and ID 130. 
In both cases the best-fit to the observed light distribution of the lens galaxy is achieved using 2 
components, i.e. an inner (more compact) Sersic profile and an exponential disk profile. These 
components are shown as a function of the distance from the galaxy center for ID81 (left-hand 
panel) and ID 130 (right-hand panel). 
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Figure S 9: Lens modelling of ID81 and ID130. The LENSMODEL software was used to fit the 
position of the peaks in the SMA images. A SIE model was assumed for the mass distribution 
of the foreground lenses. The image positions used in the fit are indicated by the red dots and 
correspond to the peaks in the SMA images (red contours in steps of 6a, 8cr, 10a etc.). The 
blue stars are the best-fit positions of the lensed background source, assumed to be point-like. 
The best-fit position of the lens galaxy and of the background source are marked by the yellow 
and the blue dots, respectively. The caustic lines and critical lines of the best-fit lens model are 
indicated in green and yellow, respectively, while the yellow dashed line shows the major axis 
of the mass model. The simulated image shows the lensed image of a background source (after 
convolution with the SMA point spread function and added noise) described by a Gaussian 
profile with FWHM= 0.2". 
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Tables 



Table S 1: Submillimeter/millimeter fluxes for the lens candidates. The official H- ATLAS 
ID, according to lAU (International Astronomical Union) convention, is that derived from the 
position of the sources in the H- ATLAS SDP catalogue. The quoted errors on the Herschel flux 
densities include statistical errors, confusion noise and an absolute calibration uncertainty of 10 
per cent for PACS/100//m, 20 per cent for PACS/160//m and 15 per cent for SPIRE. 3a upper 
limits at PACS wavelengths are provided when no detection is achieved at that significance level. 
Note that ID130 lies just outside the region covered by PACS. Fluxes at 880 /im are from the 
SMA, while those at 1200 //m are from the MAMBO array at the Institut de Radioastronomie 
Millimetrique (IRAM) 30 m telescope. 



H-ATLAS ID 


SDP ID 


Sioo 


•Sieo 


'S'250 


S350 


S500 


Ssso 


S1200 






(mJy) 


(mJy) 


(mJy) 


(mJy) 


(mJy) 


(mJy) 


(mJy) 


H-ATLAS J090740.0-004200 


9 


187±57 


416±94 


485±73 


323±49 


175±28 




7.6±0.8 


H- ATLAS J091043. 1-000321 


11 


198±55 


397±90 


442±67 


363±55 


238±37 




12.2±1.2 


H-ATLAS J090302.9-014127 


17 


78±55 


182±56 


328±50 


308±47 


220±34 




15.3±1.3 


H- ATLAS J09031 1.6+003906 


81 


<62 


<83 


129±20 


182±28 


166it27 


76.4±3.8 


20.0±0.7 


H-ATLAS J091305.0-005343 


130 






105±17 


128±20 


108±18 


39.3±2.3 


11.2±1.2 



41 



Table S 2: UV/optical/NIR photometry for the lens candidates. UV data are from GALEX, 
optical photometry is from SDSS and NIR data are from UKIDSS [as reprocessed by GAMA; 
(11)]. Where photometric measurements are not listed it means that the source is not covered at 
those wavelengths. The 3-a upper limits shown within parenthesis for the UKIDSS wavelengths 
are obtained from the residual image after the best-fit GALFIT model of the source has been 
subtracted off. These limits are used to constrain the SED of the background source. 



SDP ID 


9 


11 


17 


81 


130 


GALEX FUV (^Jy) 








0.23±0.18 




GALEX NUV (^Jy) 








1.9±1.1 




SDSS u (^Jy) 


0.24±0.23 


0.57±0.59 


3.3±1.6 


3.9±2.0 


1.7±1.7 


SDSS g (p]y) 


L79±0.43 


1.01±0.45 


3.9±6.4 


24.9±1.1 


19.41±0.72 


SDSS 1- (^Jy) 


5.81±0.70 


3.94±0.65 


7.7±1.0 


114.8±2.1 


66.1±1.2 


SDSS i (fiJy) 


14.9±L1 


11.3±1.0 


15.3±1.5 


197.7±3.6 


108.6±2.0 


SDSS z (pjy) 


27.0±3.7 


21.5±4.2 


11.8±6.0 


278.0±3.6 


143.2±6.6 


UKIDSS Y (^Jy) 






27.7±9.5(<6.6) 


321.3±3.2(<6.3) 




UKIDSS J (AiJy) 




102.4±9.8(<16) 


56±17(<12) 


367±11(<9.2) 




UKIDSS H (^Jy) 


73±15«5.0) 


237±17«14) 


107±19(<8.2) 


508.1±5.3(<8.5) 




UKIDSS K (^Jy) 


132±24(<6.5) 




108±23(<9.7) 


573.7±6.2(<14) 





Table S 3: Technical information on the SMA follow-up observations. This includes: the date 
the measurements were taken (Date), the configuration of the antennas (Conf.; VEX='very- 
extended', SUB=' sub-compact', COM='compact'), the number of antennas used (Nant.), the 
projected baselines lengths (min/mean/max Pr Baselines), the Local Oscillator Frequency (LO 
Freq.), and the on— source integration time (Int. time). 



SDP ID 


Date 


Conf. 


Nant. 


min/mean/max 
Pr Baselines (m) 


LO Freq. 
(GHz)l 


Int. time 
(min) 


81 


25FeblO 


VEX 


7 


69/281/509 


340.7 


289 


130 


28FeblO 


VEX 


7 


76/289/509 


340.7 


298 


81 


16MarlO 


SUB 


5 


6/ 17/ 25 


340.7 


144 


130 












152 


81 


09 Apr 10 


COM 


6 


9/ 38/ 69 


341.6 


153 


130 












144 



"Total bandwidth coverage is LO-8 to LO-4 (LSB) and LO+4 to LO+8 (USB) for a total of 8 GHz. Tlie small difference in LO 
Frequency between compact configuration observations and the subcompact and very extended observations is not important in this 
context. 
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Table S 4: GALFIT results for the five gravitational lens systems. The 3— cr upper limits given 
are for an extended source and derived from the distribution of 1.5 arcsecond radius aperture 
photometry of the Keck maps. 



H- ATLAS ID 


profile,^ 


•2h 


radius^ 


6i 


Axis ratio 


Anglef, 


g-3o- 


i-3a 








(arcsec) 






deg. 


(/iJy) 


(/iJy) 


9 


Sersic 


1.07 


0.85 


5.36 


0.72 


56.76 


0.162 


0.641 


11 


Sersic 


1.03 


1.10 


2.97 


0.65 


39.61 


0.229 


0.442 


17 


Sersic 


1.07 


0.61 


0.54 


0.71 


63.25 


0.202 


0.404 




Sersic 




1.36 


0.91 


0.69 


12.83 






81 


Sersic 


1.13 


0.70 


2.82 


0.78 


36.45 


0.130 


0.202 




ExpDisk 




1.20 




0.72 


0.62 






130 


Sersic 


1.00 


0.32 


1.23 


0.52 


56.82 


0.198 


0.351 




ExpDisk 




1.11 




0.55 


-54.64 







"ExpDisk = exponential disk profile 
^Reduced 

' Radius for Sersic and disk scale length for ExpDisk 
''Sersic index 

''Angle measured east of north 
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